clear; clc; close all;

% 固定参数
q = 0.4;
b = 1.2; P = 1; Q = 1;
N = 2000;        % 轨迹长度
m = 20;          % Wolf算法临近点时间间隔

% a参数范围
a_list = linspace(2.2, 2.7, 100);

% 两组初始值
X0_list = [-0.5, 0.5; 
           3,   -3];

colors = ['b', 'r'];

MLE_all = zeros(2, length(a_list));

parfor idx = 1:length(a_list)
    a = a_list(idx);
    for xi = 1:2
        x0 = X0_list(xi,1);
        y0 = X0_list(xi,2);
        [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N);
        X = [x; y];
        MLE_all(xi, idx) = MLE_Wolf(X, m);
    end
end

% 绘图
figure; hold on;
plot(a_list, MLE_all(1,:), 'b-', 'LineWidth', 1.5);
plot(a_list, MLE_all(2,:), 'r-', 'LineWidth', 1.5);
xlabel('参数 a');
ylabel('最大李雅普诺夫指数 (MLE)');
title('Fig.7(a) MLE谱，分数阶 q=0.4');
legend('x_0=(-0.5,0.5)', 'x_1=(3,-3)', 'Location', 'northwest');
grid on;
set(gca, 'FontSize', 14);

% ------------ 轨迹计算函数 -------------

% ------------ Wolf方法计算MLE -------------
function MLE = MLE_Wolf(X, m)
    N = size(X, 2);
    MLE_sum = 0;
    count = 0;

    for i = 2:N-m
        dists = sqrt(sum((X(:,i) - X).^2, 1));
        dists(i) = inf;
        idx_min = max(1, i-m);
        idx_max = min(N, i+m);
        dists(idx_min:idx_max) = inf;

        [dist_min, k] = min(dists);

        if dist_min < 1e-12 || isinf(dist_min)
            continue;
        end

        if (i+m <= N) && (k+m <= N)
            dist1 = norm(X(:,i) - X(:,k));
            dist2 = norm(X(:,i+m) - X(:,k+m));
            if dist1 < 1e-12 || dist2 < 1e-12
                continue;
            end
            Ld = log(dist2 / dist1);
            MLE_sum = MLE_sum + Ld;
            count = count + 1;
        end
    end

    if count == 0
        MLE = NaN;
    else
        MLE = MLE_sum / (count * m);
    end
end
